Binary Black Holes in Quasi-Stationary Circular Orbits 
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We propose a method of determining solutions to tlie constraint equations of General Relativity 
approximately describing binary black holes in quasi-stationary circular orbits. Black holes with 
arbitrary linear momenta are constructed in the manner suggested by Brandt and Briigmann. The 
quasi-stationary circular orbits are determined by local minima in the ADM mass in a manner 
similar to Baumgarte and Cook; however, rather than fixing the area of the apparent horizon, we 
fix the value of the bare masses of the holes. We numerically generate an evolutionary sequence of 
quasi-stationary circular orbits up to and including the innermost stable circular orbit. We compare 
our results with post-Newtonian expectations as well as the results of Cook and Baumgarte. We also 
generate additional numerical results describing the dynamics of the geometry due to the emission 
of gravitational radiation. 

PACS numbers: 95.30.Sf, 04.70.Bw, 04.20.Ex, 04.25.Dm 



I. INTRODUCTION 



A solution to the fully relativistic two-body problem 
in General Relativity is the theoretical holy grail of Ein- 
stein's theory. The advent of gravitational wave obser- 
vatories, such as LIGO have made the need for a 
two-body solution even more pressing: without accurate 
estimates of orbital parameters, investigators have little 
hope in constructing realistic waveform templates to fil- 
ter through the avalanche of observational data |^ . 

In this paper we present a method of generating so- 
lutions to the constraint equations which approximate 
a binary black hole system in quasi-stationary circular 
orbits. The quasi-stationary approximation is based on 
the existence of a helical Killing vector, which implies the 
system contains equal amounts of inbound and outbound 
gravitational radiation |^ . 

We take advantage of the conformally flat conjec- 
ture and maximal slicing, which decouples the constraint 
equations. We also adopt the "puncture" method of 
Brandt and Briigmann Q , in which our binary system is 
described by a three-sheeted Brill topology. The punc- 
ture method greatly simplifies the analysis of the sys- 
tem, discarding the need for boundary conditions at the 
throat of the holes. A numerical solution to the Hamil- 
tonian constraint is arrived at via a nonlinear adaptive 
multigrid algorithm which allows us to impose boundary 
conditions far from the punctures, while maintaining a 
relatively high density of grid points in the vicinity of 
the holes. The quasi-stationary circular orbits are deter- 
mined via a method based on a formal variational prin- 
ciple which is similar to the "effective potential" method 
of Baumgarte Q and Cook j^. The variational princi- 
ple dictates we must hold the bare masses of the holes 
fixed, and not the apparent horizon area, as we generate 
our approximately quasi-stationary configuration. This 



greatly simplifies numerical implementation, and avoids 
the need for horizon-finding algorithms. 

The paper is organized as follows: Section II is ded- 
icated to the theoretical background pertaining to the 
initial value problem, including the form of the con- 
straint equations, the puncture method, and the vari- 
ational principle we employ to generate a sequence of 
quasi-stationary circular orbits. Section III briefly de- 



scribes the numerical methods employed, and the tests 
the code was subjected to. Section ^ presents the nu- 
merical results of the code for the binary system in addi- 
tion to comparing the data to Post-Newtonian expecta- 
tions and to the results of Baumgarte and Cook. Finally, 
Section ^ gives a brief summary. 



II. THEORETICAL BACKGROUND 

A. The Initial Value Problem 

The vacuum constraint equations in the 3+1 formalism 
[Q are the cornerstone of the initial value problem. The 
Hamiltonian constraint is given by 



(1) 



where ICab is the extrinsic curvature associated with the 
three- metric 7ah, and R is the associated Ricci scalar. 
The momentum constraint is given by 



(2) 
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where _D° is the covariant derivative compatible with the 
three-metric jab- These constraint equations must be 
satisfied on every spacelike hypersurface, but do not de- 
termine the dynamics of the geometry. 

A solution to the initial value problem consists of a 
three-metric jab and an extrinsic curvature ICab which 
satisfy Eqs. and (^. However, one notes the equa- 
tions are coupled, nonlinear equations — this greatly com- 
plicates any efforts to generate either analytic or numer- 
ical solutions to the initial value equations. To simplify 
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the analysis, we assume the three-metric 7ab is confor- 
mally flat. This introduces the conformal factor tp via 
the relationship 



lab 



(3) 



where Qab is the flat three-metric. Following Bowen and 
York we have the extrinsic curvature transform ac- 
cording to 



ICab = -0 Kab- 



(4) 



We call Kab the conformal extrinsic curvature. 

In addition to the conformally flat assumption, we 
choose the maximal slicing gauge Q which corresponds 
to the conformal extrinsic curvature being trace-free: 



= 0. 



(5) 



By assuming a conformally flat three-metric, as well 
as choosing the maximal slicing gauge, we effectively 
decouple the constraint equations, making them more 
amenable to study. In particular, the Hamiltonian con- 
straint equation becomes 
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(6) 



where the operator is the flat-space Laplacian. Like- 
wise, the momentum constraint simplifies to 



(7) 



Hence, determination of Kab and ip via Eq. and 
Eq. completely specify the geometry on an initial 
space-like hypersurface. In practice, one would take this 
initial data and evolve it in time via the dynamical Ein- 
stein equations. 



B. The Puncture Method 



where is the linear momentum of the black hole mea- 
sured on the upper sheet of the geometry. Because the 
conformal metric is flat, we may use normal Cartesian 
coordinates, where n° is a radial normal vector of the 
form n° = x°/r, and = x'^ + y'^ + z'^ . For a geometry 
which consists of N black holes, the conformal extrinsic 
curvature is given by 



N 
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ah 



(9) 



i=l 



with each term in the sum having its own coordinate ori- 
gin. Note for N black holes on the upper sheet, there are 
a total of iV -I- 1 sheets in the geometry, each with its own 
separate asymptotically flat region. 

Brandt and Briigmann then proceed to solve the 
Hamiltonian constraint, given by Eq. (^), rewriting the 
conformal factor ip as 



a 



where 1 /a is given by 
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The quantity is the Newtonian mass of the z'^ black 
hole, and r(i) is the hole's corresponding location. 

With these identifications, the Hamiltonian constraint 
becomes a nonlinear elliptic equation for u: 



(12) 



where /3 is related to the conformal extrinsic curvature 
K"'' via 
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f3 = -a''K''''Kab 



(13) 



Brandt and Briigmann [Q propose a method which 
avoids some of the complications introduced by the Mis- 
ner data set, most notably the requirement of boundary 
conditions at the throat of the holes They choose a 
topology in accordance with the three-sheeted Brill data, 
and propose a method to compactify the asymptotically 
flat regions on the lower sheets. This compactification 
effectively simplifies the domain of integration in which 
the constraint equations must be solved. Simply stated, 
with the compactification of an asymptotically fiat region 
to a single point in R^, not only is the need to impose 
boundary conditions on the throat eliminated, but it also 
eliminates the need to excise domains of integration. 

Bowen and York note solutions to the momentum 
constraint, given by Eq. (|^), are well-known for a single 
black hole with arbitrary linear momentum. The trace- 
free extrinsic curvature which satisfies the momentum 
constraint is of the form 



T^ab 
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Recall we demand asymptotic fiatness on the upper 
sheet, which determines the boundary condition on u at 
spatial infinity. This Robin boundary condition B has 



du 
dr 



l-u 



(14) 



which demands that u consist of a monopole term as one 
recedes far from the black holes. 

Hence, Eq. (|) coupled to Eq. (|l|) and Eq. (|l|) com- 
plete our mathematical description of the initial value 
problem. 

In the framework of the puncture method, it is rela- 
tively straightforward to show that the ADM mass 
of a single hole as measured from that hole's asymptoti- 
cally flat region, called the hare mass to, is related to the 
Newtonian mass M of the hole measured on the upper 
sheet: 



M, 



(15) 
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where is the value of the conformal factor correction 
at the location of the i**^ puncture and is the coordi- 
nate separation distance between the i and j holes. The 
total angular momentum J and the puncture separation 
distance Dij determine each hole's linear momentum P 
via J — PDij. The values of the linear momentum P 
and the coordinate separation Dij determine the form of 
the conformal extrinsic curvature Kat via Eq. (ph. 



Once the ADM mass for the system has been calcu- 
lated for a range of separation distances D — all the while 
holding the angular momentum J and the bare masses m 
fixed — we determine the approximate location of a circu- 
lar orbit via 
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dD 



(17) 



C. The Variational Principle 

Baker and Detweiler ITlll formulate a variational prin- 
ciple for binary neutron stars which may be modified for 
binary black holes in quasi-stationary circular orbits. For 
fixed values of the bare mass m and total angular mo- 
mentum J, solutions to the constraint equations given 
by Eqs. (|^) and (|l^) generate an extremum (or station- 
ary point) in the ADM mass only when the geometry is 
a solution to the quasi-stationary Einstein equations: 



The variational principle also allows us to calculate the 
orbital angular frequency fl, given by 



SE 



ADM 



= ~2NSt 



(16) 



In the above equation, -Badm is the ADM mass of the 
system as measured on the upper sheet, N is the ratio of 
the lapse functions on the lower sheets to that of the up- 
per sheet, m is the bare mass of a hole as measured at a 
lower sheet's asymptotic infinity, is the orbital angular 
frequency of the binary system, and J is the total angular 
momentum of the binary system. Equation (|l^) assumes 
the phases of the inbound and outbound gravitational 
radiation are held fixed as one generates an evolutionary 
sequence of quasi-stationary circular orbits. 

The application of the variational principle for binary 
black holes follows. Specify the value of the total angular 
momentum J and the bare masses m. These quantities 
are to be held fixed during the analysis. Then, initialize 
the puncture separation distance D. 

For some initial guess of the Newtonian mass M, solve 
the Hamiltonian constraint equation given by Eq. dl^ ) 
subject to the Robin boundary condition given by 
Eq. (H). Once u is known, use it to determine new val- 



ues for the Newtonian masses, as determined by Eq. (15), 
so as to hold rn fixed. 

Repeat the above procedure by re-solving the Hamil- 
tonian constraint equation using the newly determined 
Newtonian masses, iterating until the Newtonian masses 
have converged to within some predetermined tolerance 
and give the appropriate value of the bare mass. Only 
then assume that the Hamiltonian constraint is solved 
subject to fixed values of the total angular momentum J 
and the bare mass m. At this point, calculate the ADM 
mass for the system. 

Next, decrement the hole separation distance D, hold- 
ing the total angular momentum J and the bare masses 
m fixed to their initial values. Repeat this procedure for 
a range of hole separations, calculating the ADM mass 
at each separation distance after reasonable convergence 
in the Newtonian masses has been obtained. 



dEADM 



dJ 



(18) 



In fact, the variational principle implies that if the ini- 
tial data are within some small measure S of true quasi- 
stationary data, then E'adMj Ji and N are within 
of their true quasi-stationary values. 

Note that our procedure is similar to those of Baum- 
gartc and Cook but in their heuristic effective 
potential method they hold the area of the apparent hori- 
zons fixed as they generate a member of their evolution- 
ary sequence of circular orbits. 



III. NUMERICAL METHODS 

We briefly discuss the numerical methods imple- 
mented to solve the nonlinear constraint equation, given 
by Eq. (p^). We developed an algorithm in C-f-|- to 
solve the Hamiltonian constraint via adaptive multigrid 
methods. We chose the multigrid method because of its 
relative ease of implementation, as well as its speedy 
convergence rate [|l^, |l3|. The adaptive portion of 
the algorithm allows the Robin boundary condition 
to be imposed arbitrarily far from the holes, subject 
only to machine memory limitations. This is achieved 
by starting with a cubical grid comprised of grid 
points, with a physical grid length of L. Each successive 
adaptive grid is also comprised of it' grid points, but 
has a physical grid length twice that of the previous 
adaptive grid. This results in a series of embedded 
grids that yields a high grid density in the area around 
the holes, while allowing the boundary condition to be 
imposed far from the holes. Coupled to the symmetry 
of the equal mass scenario, this allowed for a relatively 
large number of total grid points to be employed. For 
instance, if the smallest grid has n grid points on a side 
and a total of N adaptive levels, the total number of 
grid points employed in the calculation is on the order 
of Nn^. Detailed discussion of the adaptive multigrid 
algorithm, as well as the complete C++ computer code, 
may be found in p4|. 
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FIG. 1: The radial dependence of the normahzed energy den- 
sity r^p(r)l(PIMf versus the normahzed distance r/M. The 
quantity p(r) is defined as l/27r J /3(1 + au)'"^ dO.. 



Numerical Tests 



As demonstrated in [Q, the ADM mass may be deter- 
mined via 



DM 



1 

+ 2^ 



a, 



-7j3 



(19) 



It is straightforward to analyticaUy determine the 
ADM mass for a single boosted black hole in the limit 
of small P/M. Restricting the analysis to u ^ 1 and 
P/M < 1.0, we find the result to be 



DM 



M H . 

8 M 



(20) 



One's attention is immediately drawn to the factor of 
5/8 in front of the "kinetic" term. One would prefer the 
initial value solution described by Eqs. (|^)-(|l^) would 
yield the expected Newtonian factor of 1/2. That it does 
not reflects the fact that even for this simple example, 
the initial data appears to contain mass-energy in ad- 
dition to that of the black hole itself. Figure |^ shows 
the radial dependence of the energy density; integration 
of this curve from r = to r = oo yields the factor of 
5/8. It is impossible to pinpoint the location of this ad- 
ditional energy; all that can be said is the initial data 
contains additional energy "frozen" into the geometry. 
Nonetheless, we may use this geometry as a test bed of 
the code, and we will ultimately apply the variational 
principle to the geometry to examine binary black holes 
in quasi-stationary circular orbits. 

In numerical runs of the single boosted hole with 
P/M = 0.1 and a physical grid length for the small- 
est grid of 4M, we found agreement with the theoretical 



ADM mass correction to be 1.3%, 0.6%, and 0.4% for cu- 
bical grids with n = 11, n = 19, and n = 35 grid points 
on the side of the smallest adaptive grid, respectively. 
These grid resolutions translate to 2, 4, and 8 grid points 
across the diameter of the throat of the hole, respectively. 
For all three resolutions, the Robin boundary condition 
was imposed at a distance of approximately 4000M from 
the puncture. 

A second test of the code, as described in Q, is two 
boosted black holes located on an axis at coordinate dis- 
tances ±1? and moving towards each other. The theoret- 
ical expression for the ADM mass is not as straightfor- 
ward to derive as for a single hole. Once again, assume 
the correction to the conformal factor w « 1, and limit 
the analysis to P/M <?C 1.0. After Taylor expanding the 
second integrand of Eqn. ( p^ in factors of D/M, we find 
the ADM mass for the system to be 



Ea 



DM 
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35 


\m 



(21) 

Note Eqn. ( |2l] ) is accurate only up to fourth order in 
the expansion parameter D/M. In numerical runs with 
P/M 0.1, D/M = ±0.1, and a physical grid length 
for the smallest grid of 4M, we found agreement with 
the theoretical ADM mass correction to be 5.4%, 1.7%, 
and 0.9% for cubical grids with n = 11, n = 19, and 
n = 35 grid points on the side of the smallest adaptive 
grid, respectively. Once again, for all three resolutions, 
the Robin boundary condition was imposed at a distance 
of approximately 4000Ajf from the punctures. 



IV. NUMERICAL RESULTS 

We now present the numerical results from our adap- 
tive multigrid computer code used to model binary black 
holes. Only binary systems with black holes having zero 
intrinsic spin were investigated. 

Figure || is an example of the curv es generated by the 
procedure described in Section II C , for fixed values of 
the total angular momentum J and bare mass m. The 
vertical axis is the ADM mass of the system, and the hor- 
izontal axis is the coordinate separation distance between 
the two holes. The location of the circular orbits is de- 
noted by the dashed curve, with diamonds representing 
the locations of the circular orbits for particular values 
of J. This particular figure was generated for a grid with 
n = 11 grid points on the side of the smallest adap- 
tive grid, with the Robin boundary condition imposed at 
approximately 10001?, where D is the coordinate sepa- 
ration distance of the holes. The bottommost solid line 
curve corresponds to a value of the angular momentum 
J ~ 2.93m^, and the topmost solid line curve corresponds 
to a value of J = 3.00rn^. The evolutionary sequence 
of circular orbits terminates at a coordinate separation 
distance oi D = 5.8395m and total angular momentum 
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FIG. 2: The ADM mass curves for fixed values of J and 
m, and the corresponding evolutionary sequence of circular 
orbits for n — 11 grid points on a side of the smallest adaptive 
grid. The Robin condition is imposed at a distance of lOOOD. 
The values of J range from 2.93m^ for the bottommost solid 
line curve to S.OOm'^ for the topmost solid line curve. Our 
approximation of the evolutionary sequence terminates at J = 
2.937m^, the innermost stable circular orbit. 



J = 2.937m'^, which corresponds to the innermost stable 
circular orbit. 

The resulting evolutionary sequence of circular or- 
bits, denoted by the dashed line, can be viewed in the 
following fashion. A realistic binary system will radiate 
gravitational waves, causing their orbital separation dis- 
tance to decrease slowly. The binary system will progress 
along the evolutionary sequence, until the binary system 
reaches the innermost stable circular orbit. It is after 
the innermost stable circular orbit that we may no longer 
employ our variational principle, simply due to the lack 
of stable circular orbits for smaller separation distances. 
Because of this, we can not determine the dynamics of 
the final inward spiral and plunge. Despite this short- 
coming, the variational principle, and the evolutionary 
sequence it generates, yields valuable information about 
the evolution of the system up to and including the in- 
nermost stable circular orbit. 

One problem that arises in the analysis is the question 
of black hole mass. Specifically, what does one mean 
when one speaks of the mass of a black hole in a binary 
system? There is no general conse nsus on the answer to 
this question. Recall from Section II B that we encoun- 



tered two types of mass when describing the punctures of 
Brandt and Briigmann: the bare mass m and the New- 
tonian mass M . 

However, these masses are not the only used to de- 
scribe black holes. For instance. Cook |^ and Baumgarte 
1^ choose to use the mass associated with the area of the 
apparent horizon, described by Christodoulou's formula 

£1- 

A different measure of the mass of a black hole is the 



rest mass Q , described by the well-known formula from 
special relativity: 



^ADM 



(22) 



where E'adm is the ADM mass of an isolated black hole 
with linear momentum P. We have opted to use the rest 
mass when comparing our numerical results to the ex- 
pected post-Newtonian results, simply because its value 
depends only upon the properties of the single hole, and 
not the interaction energy between the holes. However, 
Ej-est includes the additional energy which is "frozen" into 
the initial data, and would be expected to effect the dy- 
namics of the system. For a binary system with each hole 
having a linear momentum P, the rest mass of each hole 
is calculated by numerically evaluating the ADM mass 
for a single, isolated, boosted black hole with linear mo- 
mentum P. 

Kidder et al. provide (post)^-Newtonian expres- 
sions for the total angular momentum J, orbital angular 
frequency Q, and the binding energy E\^. Instead of cal- 
culating the binding energy based on the mass associated 
with the apparent horizon ^, we define the binding 
energy to be 



Eh = Eadm — 2£^r( 



(23) 



We also introduce the reduced mass /i and the total mass 
m, defined via 



and 



(24) 



(25) 



The (post)^-Newtonian equations, written in the form 
of Cook P, are 



Eb 



Eb 



1 / fj,m\'^ 

2 V~/ 



1 + 
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(26) 



1069 
384 



(mr!)4/3 



(27) 



and finally 



(—] ^ (mr!)-2/3 



+ ^(ml7)^/3 + 



.(28) 



We now present the figures of our numerical results, 
compared to (post)^-Newtonian expectations. 
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FIG. 3; The normalized binding energy i?b/M versus the nor- 
malized angular momentum J/fim. 
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FIG. 4: The normalized binding energy Eh/i-t versus the nor- 
malized angular frequency mfi. 



Figure ^ displays the normalized binding energy / ^ 
versus the normalized angular momentum J/fj.m. The 
solid line is the theoretical prediction given by Eq. (26). 
One notes that the curves corresponding to n = 19 and 
n = 35 actually cross the theoretical prediction, as op- 
posed to asymptotically approaching the solid line. We 
believe that if the grid resolution were to be increased, 
the resulting data would conform more closely to the 
post-Newtonian prediction for large J/ fim. 

Figure ^ displays the normalized binding energy 
i^b/A* versus the normalized angular frequency mfi. The 
solid hne is the theoretical prediction given by Eq. (27). 
Note the Newtonian regime of the numerical curves, 
corresponding to smaller values of mQ, agree with the 
post-Newtonian prediction. There is agreement between 
the theoretical prediction and the two higher resolution 
curves up until mfl « 0.08. At this point, the curvature 
of the two higher resolution curves deviates from the the- 
oretical prediction, and tends to more negative values of 
the normalized binding energy. 

Figure ^ displays the normalized angular momentum 
J/ lim versus the normalized angular frequency mfl. The 
theoretical prediction is given by Eq. (p8|), and is denoted 
in the figure by the solid line. The curves indicate con- 
vergence at both the Newtonian and the relativistic ends 
of the figure. 

Before we compare the data to that of Baumgarte ||] 
and Cook we feel it is important once again to stress 
that there is some ambiguity as to what one means when 
one talks about the mass of a black hole in a binary sys- 
tem. Both Cook and Baumgarte associate the area of the 
apparent horizon to the mass of an individual hole, and 
hold the area of the apparent horizon fixed as they gener- 
ate their evolutionary sequence of circular orbits. We, on 
the other hand, choose to use the rest mass to describe 
the individual hole, and held the bare mass of each hole 
fixed as we generated our evolutionary sequence of circu- 



Post-Newtonian 
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FIG. 5: The normalized angular momentum J/fj,m. versus the 
normalized angular frequency mQ. 



lar orbits. 

We begin by noting that we cannot compare our coor- 
dinate separation distance D /m to the proper separation 
distances determined by Baumgarte and Cook. Our main 
limitation in this is our lack of knowledge concerning the 
locations of the apparent horizons of the holes. How- 
ever, we may still get an approximate measure d of the 
separation distance of the holes at the innermost stable 
circular orbit: we assume the quadrupole moment of the 
system, as measured at infinity, consists of two objects 
of mass ^-Eadm separated by the distance d. We present 
the innermost stable circular orbit's normalized coordi- 
nate separation D/m, as well as the separation distance 
d normalized by the system's ADM mass, in Table | 

Table H compares the results for the innermost sta- 
ble circular orbit for the highest resolution of this work 
(n = 35 grid points on the side of the smallest adaptive 
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TABLE I: The normalized coordinate separation distance 
D/m, and the normalized separation distance (I/Eadm for 
three resolutions at the innermost stable circular orbit. 



Resolution 


D/m 


d/ -Badm 


n = 11 


2.915 


3.078 


n = 19 


2.143 


2.290 


n = 35 


2.195 


2.349 



TABLE II: Comparison of the innermost stable circular orbit 
data for Cook, Baumgarte, and this work. 



n=35 

Newtonian 




Data 


£b//i 


J / lira 


mO. 


Cook 


-0.09030 


2.976 


0.172 


Baumgarte 


-0.092 


2.95 


0.18 


This work 


-0.10792 


2.937 


0.179 



grid) to that of Cook |^ and Baumgarte [||, as listed in 
Baumgarte. 

The agreement between our data and that of Cook 
and Baumgarte may lend some credence both to their 
method of approach, as well to the variational principle 
we employed. On one hand, our method is based on a 
mathematically derived variational principle for binary 
black holes in quasi-stationary circular orbits in which 
the variational principle dictates we must hold the bare 
mass fixed, and not the area of the apparent horizon. 
This greatly simplifies the analysis as one does not need 
to be concerned with the location of the apparent hori- 
zon, which can be a daunting numerical task ||l7| , 
On the other hand, the agreement of Cook and Baum- 
garte's results with ours may indicate that there is some 
interesting physics yet to be discovered concerning the re- 
lationship between the variational principle we employed 
and the methods employed by Cook and Baumgarte. 



A. Gravitational Waves 

An algorithm was developed to calculate the 
quadrupole moment of the field, which yields the di- 
mensionless luminosity of the gravitational waves via the 
quadrupole moment formula |l9[, and subsequently de- 
termines the dynamics of our binary system. 

Figure ^ displays the gravitational wave luminosity L 
as a function of the normalized separation distance D/m. 
The three different resolutions are displayed, as well as 
the Newtonian expectation for the gravitational wave lu- 
minosity. The numerical results appear to be converging 
to a limit which is in line with the Newtonian expecta- 



FIG. 6: The dimensionless gravitational wave luminosity L as 
a function of the normalized separation distance D/m. 



tion for large separation distances. Also, the numerical 
results are in excellent agreement with each other as the 
separation distance approaches the innermost stable cir- 
cular orbit. 

The quasi-stationary approximation is a powerful tool, 
but it is natural to inquire about its range of validity. 
Specifically, can the information extracted from the grav- 
itational radiation indicate a point in which we may no 
longer rely upon the quasi-stationary approximation to 
describe our binary system? 

We begin by noting the gravitational wave luminosity 
L may be used to determine the time derivative of the 
orbital angular frequency via 



dn _ f dE ADM \ ^ ^ 

dt ~ \ dn J 



(29) 



The results of this calculation for the three resolutions 
under study are shown in Fig. |^ which shows m?dVl/dt 
versus the normalized orbital angular frequency mf2. As 
indicated in the figure, as the system approaches the in- 
nermost stable circular orbit, the rate of change of the 
angular frequency diverges. 

The information from dQ/dt may in turn be used to 
yield the orbital angular frequency as a function of time. 
This is achieved via 



t — hsco 



dn 

f^ISCO \ dt ) 



(30) 



The quantity t — tisco is a negative number, and is a 
measure of the amount of time that must elapse before 
the system reaches the innermost stable circular orbit. 
Figure || shows the normalized orbital angular frequency 
niQ versus the normalized time scale {t — tisco) /'ti. The 
solid bold curve is the (post)^-Newtonian equation of 
Blanchet et al. The diamond on each numerical 

curve corresponds to the time when the system begins 
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TABLE III: Comparisons between the quasi-stationary ISCO 
and the beginning of the last orbit, or BLO. 



FIG. 7: The normahzed rate of change of the orbital angu- 
lar frequency m^dfl/dt versus the normalized orbital angular 
frequency niQ. As the system approaches the innermost sta- 
ble circular orbit, the orbit evolves so rapidly that m?dQ,/dt 
effectively diverges, as indicated by the nearly vertical asymp- 
totes. 



n-19 
- n-35 
. 2PN 




FIG. 8: The normalized orbital angular frequency mf2 versus 
the normalized time scale — tisco) /m. For a particular 
data set, the diamond on each curve marks tsLO, the approx- 
imate time for one more complete orbit before reaching the 
innermost stable circular orbit. 



its last orbit, ieLOi before reaching the innermost stable 
circular orbit; i.e., 



^I^BLO — ^ISCol 

2^r 



(31) 



is satisfied. 

In Table III, we list various parameters associated 
with the beginning of the last orbit, or "BLO", as deter- 
mined by Eq. (^ij). We compare this with the values of 
the innermost stable circular orbit, or "ISCO", as deter- 
mined by the evolutionary sequence of quasi-stationary 





rnQ, 




D/m 


v/c 


ISCO 


0.179 


2.937 


2.19 


0.56 


BLO 


0.049 


3.401 


6.62 


0.25 



circular orbits. All values are taken from the highest 
resolution data. Table III shows the parameters of the 



BLO are far removed from the parameters associated 
with the ISCO. However, the value oi v/c « 0.25 for 
the BLO indicates the quasi-stationary method probes 
regimes where post-Newtonian approximations may start 
to break down. 



V. SUMMARY 

In this paper, we propose a method of generating se- 
quences of approximate quasi-stationary circular orbits 
for binary black hole systems similar to the method of 
Cook Q and Baumgarte Q . The method is based on the 
variational principle of Baker and Detweiler ||Tl| , and em- 
ploys the "puncture" method of Brandt and Briigmann 

i- 

We developed an adaptive multigrid algorithm which 
solved the nonlinear Hamiltonian constraint and allowed 
us to extract useful orbital parameters for the system, 
up to and including the innermost stable circular orbit. 
Despite the additional energy content of the geometry, 
our results agree with post-Newtonian expectations, and 
are in line with the work of Baumgarte and Cook [|j . 
We also determined, via extraction of gravitational wave 
information, the beginning of the last orbit, or BLO. 

Our method is straightforward as holding the bare 
mass fixed does not require locating the apparent horizon 
of the holes. Coupled to the puncture method of Brandt 
and Briigmann, numerical implementation via adaptive 
multigrid methods is relatively easy. 

This paper reveals the power of the variational princi- 
ple. Relaxation of the conformally flat conjecture is the 
logical next step. This will couple the constraint equa- 
tions and introduce complications in determination of a 
solution, but it will eliminate the additional energy con- 
tent of the geometry, and allow for more realistic simu- 
lations of quasi-stationary binary black holes. 
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